

library(tidyverse)
library(rio)
library(countrycode)
library(wesanderson)

setwd("~/06_goldstein_rivers_tomz_2007")

full_data <- import("replication_data.dta")


world <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  mutate(wdicode = countrycode(region, 'country.name', 'wb', warn = T))

data_og1 <- full_data %>%
  mutate(wdicode = countrycode(name_1, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,gdp,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique() 

data_og2 <- full_data %>%
  mutate(wdicode = countrycode(name_2, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,gdp,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique() 

data_og <- data_og1 %>%
  full_join(data_og2) %>%
  group_by(wdicode) %>%
  mutate(n = sum(n)) %>%
  unique() %>%
  right_join(world) %>%
  cbind(name = "Original Data")


data_wdi1 <- full_data %>%
  mutate(wdicode = countrycode(name_1, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant2000_wdi2005,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique()

data_wdi2 <- full_data %>%
  mutate(wdicode = countrycode(name_2, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant2000_wdi2005,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique()

data_wdi <- data_wdi1 %>%
  full_join(data_wdi2) %>%
  group_by(wdicode) %>%
  mutate(n = sum(n)) %>%
  unique() %>%
  right_join(world) %>%
  cbind(name = "WDI 2005")

data_pwt1 <- full_data %>%
  mutate(wdicode = countrycode(name_1, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant1996_pwt61,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique()

data_pwt2 <- full_data %>%
  mutate(wdicode = countrycode(name_2, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant1996_pwt61,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique() 

data_pwt <- data_pwt1 %>%
  full_join(data_pwt2) %>%
  group_by(wdicode) %>%
  mutate(n = sum(n)) %>%
  unique() %>%
  right_join(world) %>%
  cbind(name = "PWT 6.1")

data_maddison1 <- full_data %>%
  mutate(wdicode = countrycode(name_1, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant_maddison2003,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique() 

data_maddison2 <- full_data %>%
  mutate(wdicode = countrycode(name_2, 'country.name', 'wb', warn = T)) %>%
  select(wdicode,imports,gatt_1,gatt_2,ptarecip,ptanonrecip,gsp,currencyunion,colorbit,log_gdp_constant_maddison2003,distance,share_language,
         share_border,landlocked,island,land_area) %>%
  drop_na() %>%
  group_by(wdicode) %>%
  tally() %>%
  unique() 

data_maddison <- data_maddison1 %>%
  full_join(data_maddison2) %>%
  group_by(wdicode) %>%
  mutate(n = sum(n)) %>%
  unique() %>%
  right_join(world) %>%
  cbind(name = "Maddison")

map_data <- rbind(data_og,data_wdi,data_pwt,data_maddison) %>%
  mutate(name = as.factor(name),
         name = relevel(name, ref = "PWT 6.1"),
         name = relevel(name, ref = "WDI 2005"),
         name = relevel(name, ref = "Original Data"))

pal <- wes_palette("Zissou1", 100, type = "continuous")
map_data %>%
  ggplot(aes(long, lat, group = group)) + 
  geom_polygon(aes(fill = n), color = "gray92", linewidth = 0.001) + facet_wrap (~ name) + 
  theme_void() + labs(fill="Number of \nAvailable \nObservations") +
  scale_fill_gradientn(colours = pal, na.value="gray92") + 
  theme(plot.margin = unit(c(0,0,0,0), "cm")) +
  theme(strip.text.x = element_text(size = 12, face = "bold", hjust = 0.5))

# ggsave(filename = "map_available_goldstein.pdf", height = 7, width = 12)
